function P_AFAD=Pclogit_AFAD(Y_AFAD,X_AFAD,Z_AFAD,b_AFAD)

choice=2:13;
ns=numel(choice);

Y=sortrows(Y_AFAD,[1,2]);
Y=Y(:,3:end);
Y=Y-1;  %!!!### UPDATED on 06/03/2020: Y MUST be from 1 to J!!!

X=sortrows(X_AFAD,[1,2]);
X=[ones(size(X,1),1),X(:,3:end)];

Z_AFAD=sortrows(Z_AFAD,[1,2,3]);
Z_AFAD(Z_AFAD(:,3)==1,:)=[];
Z=zeros(size(Y,1),size(Z_AFAD,2),ns);
for i=4:size(Z_AFAD,2)
    Z(:,i,:)=reshape(Z_AFAD(:,i),ns,[])';
end
Z(:,1:3,:)=[];

baseAlt=1;

P_AFAD = pclogit(b_AFAD,Y,X,Z,baseAlt);








